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ABSTRACT 

Many astronomers have speculated that the solar system contains undiscov- 
ered massive planets or a distant stellar companion. The acceleration of the solar 
system barycenter can constrain the mass and position of the putative compan- 
ion. In this paper we use the most recent timing data on accurate astronomical 
clocks (millisecond pulsars, pulsars in binary systems and pulsating white dwarfs) 
to constrain this acceleration. No evidence for non-zero acceleration has been 
found; the typical sensitivity achieved by our method is a & /c ~ a fewxlCT 19 
s" 1 , comparable to the acceleration due to a Jupiter- mass planet at 200 AU. The 
acceleration method is limited by the uncertainties in the distances and by the 
timing precision for pulsars in binary systems, and by the intrinsic distribution 
of the period derivatives for millisecond pulsars. Timing data provide stronger 
constraints than residuals in the motions of comets or planets if the distance to 
the companion exceeds a few hundred AU. The acceleration method is also more 
sensitive to the presence of a distant companion (> 300 — 400 AU) than existing 
optical and infrared surveys. We outline the differences between the effects of 
the peculiar acceleration of the solar system and the background of gravitational 
waves on high-precision timing. 

Subject headings: gravitational waves — pulsars: general — solar system: general 

1. Introduction 

For over 150 years, astronomers have wondered whether the solar system contains undis- 
covered massive planets or a distant stellar companion. Small residuals in the orbit of Uranus 
motivated the search for high parallax objects near the ecliptic that resulted in the discovery 
of Pluto in 1930 (Tombaugh 1961). Numerous trans-Neptunian bodies have been discovered 
since, e.g., using the Caltech Wide Area Sky Survey (Trujillo & Brown 2001), but all these 
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objects (including Pluto) are too light to affect the major planets' orbits at the current mea- 
surement precision. Standish (1993) argued persuasively that the claimed residuals in the 
orbit of Uranus were likely due to systematic errors, such as an underestimate of Neptune's 
mass by 0.5%, rather than to an unseen massive object. The various constraints on the 
mass and position of a putative Planet X or a distant stellar companion were summarized by 
Tremaine (1990) and Hogg et al. (1991), who concluded that there was no dynamical evidence 
for an unseen companion down to the observational precision at the time of writing. 

A massive companion would exert an acceleration on the known solar system barycenter. 
Such acceleration affects the observed value of the rate of the period change of astronomical 
clocks, such as pulsars and pulsating white dwarfs. In this paper we use the most up-to-date 
timing data to look for evidence of this acceleration. 

We consider a clock far outside the known planets' orbits, at distance d in the direction 
given by unit vector n and moving with velocity v (all values are relative to the barycenter 
of the known solar system). Due to some physical process, the clock has an intrinsic period 
P which changes at a rate P. The observed period and period change rate are: 




where vt is the transverse velocity (vt = f^d, where fi is the proper motion) and a is the 
acceleration of the clock relative to the solar system barycenter. The component of this 
acceleration due to the solar system companion is denoted — a Q throughout this paper (so 
that the acceleration of the solar system due to the companion is a ). The neglected terms 
include the special relativistic time dilation which is of order v 2 /c 2 . In equation (2) the 
intrinsic period P can be replaced by the observed period P' to the same order. 

A planet of mass m p at distance d p would exert an acceleration of 
a Q _Gm p _ 18 , / m p \ ( d v \ 2 

T~^dJ~ 2xW s \mJJ\iooW) ' W 

where Mj is the mass of Jupiter; we consider timing of objects that can provide similar or 
better constraints. In Section 2 we consider constraints on the acceleration in a statistical 
sense from timing of large groups of objects. In Section 3 we consider acceleration constraints 
from individual objects. In Section 4 we briefly discuss effects of a possible background of 
low-frequency gravitational waves on timing of astronomical clocks. We summarize our 
results in Section 5. 
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2. Statistical constraints on the acceleration 

2.1. Statistical approach 

Pulsars are among the most stable astronomical clocks, their periodic signals being due 
to the rotation of a compact neutron star. About 1500 pulsars are known. Typically, the 
periods of pulsars slowly increase with time as they lose their rotational energy to particles 
and radiation. Whenever a secular period decrease is observed, it is believed that the ob- 
served value of P' is not intrinsic but rather due to the acceleration of the pulsar relative 
to the solar system. In support of this view, negative period derivatives have so far only 
been observed in globular clusters (Freire et al. 2001), arising from the pulsars' accelerations 
in the clusters' potentials. In most cases, extracting the acceleration using equation (2) is 
problematic, since (i) the intrinsic period derivative is not known (first term of eq. 2 on the 
right-hand side); (ii) line-of-sight velocities of pulsars have not been measured (second term); 
(iii) distances to pulsars are highly uncertain (Taylor & Cordes 1993; Brisken et al. 2002); 
and (iv) proper motions have been reliably measured only for relatively nearby objects. 

We now proceed to estimating typical contributions of different terms to the observed 
period change rate in equation (2). If we assume for a moment that the relative acceleration 
is at the level of a/c ~ 10~ 18 s _1 , then for ordinary radio pulsars (P ~ 1 s, P ~ 2 x 10~ 15 , 
v ~ 100 km s _1 , r ~ 1 kpc) the typical values are 

P'/P ~ 1 + 3 x 10~ 4 + 5 x 10" 4 + 5 x 10~ 4 (4) 

(in the same order as in eq. 2). The observed period change is dominated by the intrinsic 
value, and all other terms are tiny corrections of similar magnitudes. Thus ordinary radio 
pulsars are not useful for constraining the solar system acceleration on the order a /c < 10~ 18 
s _1 . An additional complication is that normal pulsars concentrate strongly toward the 
Galactic plane and therefore cannot be used to constrain the component of the acceleration 
perpendicular to the plane. 

The period distribution of radio pulsars is bimodal (Phinney & Kulkarni 1994), with a 
minimum at about 25 ms; there are only a few tens of objects with shorter periods. These 
millisecond pulsars (MSPs) typically have P ~ 10~ 2 s, P ~ 4 x 10~ 20 , v ~ 100 km s _1 , and 
r ~ 1 kpc. In this case, 

P'/P ~ 1 + 3 x 10~ 4 + 3 x 10" 1 + 2.5 x 10" 1 . (5) 

The proper motion and the possible acceleration-related corrections are comparable to the 
intrinsic value of the period derivative, whereas the radial velocity term can be neglected. 
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In Section 2.3 we implement a statistical approach to determining the solar system 
barycenter acceleration using MSPs. This approach is based on the plausible assumption 
that the intrinsic P distribution is independent of the position on the sky. Therefore, if 
the solar system accelerates relative to the MSP population, there would be a systematic 
dependence of the observed P' on the position of the pulsars in the Galaxy (this idea was 
first explored by Harrison 1977 for a small sample of ordinary pulsars). Such a change can 
be detected, provided that it is not too small compared to the dispersion of values P' in the 
population. Therefore, in contrast to ordinary pulsars, MSPs offer the hope of measuring a 
solar acceleration. 



2.2. Data handling and limitations 

From the ATNF pulsar database 1 (Manchester et al. 2005) as of January 2005, we 
selected all objects with P' < 25 ms. We rejected all objects associated with globular 
clusters (as per Harris 1996) since their accelerations relative to the solar system are likely 
to be dominated by the potentials of the systems they are residing in. We also neglected the 
young pulsar PSR J0537— 6910 from the Large Magellanic Cloud. We included only objects 
for which P' has been measured, producing a sample of 48 objects. Their distribution on the 
sky is shown in Figure 1 which also introduces the sky visualization used throughout this 
paper. 

Many of the parameters given in the database such as the period rate changes, proper 
motions, parallaxes and orbital parameters (for those MSPs that are in binary systems) are 
obtained using the analysis of the times of arrival of the pulses, as described by Kramer 
(2004). Of the 48 MSPs that we use, 35 are members of binary systems, and the period 
derivatives that we use have been corrected for the orbital motion. The relative errors of 
the period derivative measurements are typically better than 0.1%, even for MSPs in binary 
systems, so they are a negligible source of uncertainty. Only eight MSPs in our sample have 
measured parallaxes, with a typical accuracy 30%. Therefore, for most of the objects the 
distances were estimated using the measured values of the dispersion measure and the model 
of the distribution of the free electrons in the Galaxy by Taylor & Cordes (1993). These 
estimates can be uncertain by up to 40% (Brisken et al. 2002); we discuss the effects of these 
uncertainties in Section 2.3. 28 out of 48 MSPs have proper motions measured as part of 
the timing analysis. 

All values listed in the database are corrected for the effects of the major bodies of 
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the solar system and are given relative to the known solar system barycenter. Therefore, 
the accuracy of the solar system ephemerides limits the accuracy of the period derivative 
(and therefore of the relative acceleration). Fairhead (1989) and Damour & Taylor (1991) 
estimate this accuracy to be cr eph /P = 1.2 x 10~ 20 s -1 , where a ep h is the accuracy of the 
period derivative. The relative acceleration between a pulsar and the solar system that 
enters equation (2) is composed of their respective accelerations in the Galactic potential 
and their peculiar accelerations due to interactions with neighboring objects. Were there a 
companion to the solar system, the acceleration a Q it exerts must be greater than the peculiar 
acceleration of the solar system due to field objects in order to be detectable. The peculiar 
accelerations due to the field stars and giant molecular clouds in the solar neighborhood 
do not exceed (JfieujP < 1 x 10~ 20 s^ 1 (Damour & Taylor 1991). We explicitly take into 
account the relative acceleration due to the overall Galactic potential as described in Section 
2.3. 

Finally, a companion to the solar system would be displaced as a result of its orbital 
motion on the time scale that the pulsar timing exists (30+ years, or 20 years in the case 
of MSPs) leading to a change of the direction of the acceleration exerted on the solar sys- 
tem. This change is quite significant at the distance of Pluto, with angular motion of 
30deg x (dp/ '40 AU) 3 ^ 2 (t/20 years), but becomes negligible for distances to the compan- 
ion d p > 100 AU, or if we primarily use more recently discovered MSPs. Assuming one of 
these latter conditions holds, we neglect the change of the direction of the acceleration in all 
subsequent discussion. 

2.3. Millisecond Pulsars 

For each direction on the sky (a, 5) we construct two datasets based on the observational 
data for the 48 MSPs described above (in practice we used a grid spaced by 10 degrees in 
right ascension and declination, as shown in Figure 1). The first dataset contains -P//-P/ for 
all MSPs and the second dataset consists of cos#j, where 0j is the angle between (a, 5) and 
the direction to the object (ai,6i). If the solar system accelerates with a in the direction 
(ex, 5), from equation (2) we know that for all objects there is a term —a Q cos6/c which 
contributes to the observed P'/P' value. If there is no acceleration in this direction, there 
should be no correlation between the observed timing properties of MSPs and their position 
on the sky (the null hypothesis). 

The presence of the correlation is assessed using the Spearman rank correlation coef- 
ficient which is calculated after replacing all values in a dataset by their ranks (e.g., the 
minimum value is replaced by 1 and the maximum by N if the dataset consists of iV values). 



- 6- 



The advantage of using the rank sets is that the results do not depend on the distribution 
from which the datasets were drawn (Press et al. 1992). The rank correlation coefficient r s 
is a signed quantity (positive if one value on average increases as a function of the other 
and negative if it decreases); values of r s around imply there is no correlation, whereas 
values close to 1 or — 1 indicate a strong correlation. The significance of the rank correlation 
coefficient can be assessed using the probability H(r s ) that a given value of r s is obtained 
for two uncorrelated sets of data. 

For each direction we evaluate the rank correlation coefficient r s (a, 5) between the two 
datasets P[jP[ and cos#j. We then evaluate the probability of the null hypothesis U(a,5) 
(i.e., the probability that P[jP[ and cos 6^ are uncorrelated) and use these values to construct 
a "probability map" . An example of a probability map is given in Figure 2 for a simulated 
dataset. For this simulation we used the real positions of MSPs, but assigned them all the 
same P and P values and added an artificial acceleration to the data producing slightly 
different observed P' values. The probability map for the negative declinations is fully 
determined by the probability map for positive declinations, since the calculated probabilities 
by construction are invariant with respect to reversing the direction, U(a, 5) = n(7r + a;, —5). 

To produce the probability map using the real data on MSPs we corrected the observed 
period derivatives for the proper motion effect and the relative acceleration in the Galac- 
tic potential using equation (2). To this end, we constructed a new dataset consisting of 
Pi, new/Pi = P[lP[ ~ /J'idi/c — dgaij/c, where Hi are the observed proper motions of the MSPs 
(consistent with zero for 20 out of 48 objects), di are their distances from the sun and a ga i^ 
are the accelerations of the pulsars relative to the solar system in the Galactic potential pro- 
jected onto the direction to the pulsars. The typical values of the proper motion correction 
are 

%^ = ^ = 2.5x 1(T 19 s- 1 ( ^ r V ( -4-] ■ (6) 

P c \10 mas year i / \1 kpc/ 

The correction for the relative Galactic acceleration is performed using the Galactic potential 

in the form suggested by Paczynski (1990) (hereafter 'P90 model') which consists of a disk, 

a spheroid (bulge) and a spherically symmetric halo; the model is normalized to have the 

rotation velocity of 220 km s _1 at 8 kpc. This model of the Galactic potential was updated 

by Dehnen & Binney (1998). Numerically, the P90 and the Dehnen & Binney (1998) models 

produce similar results (Sun & Han 2004), but the P90 potential is easier to use since it 

has an analytic form. The typical values of this correction are a few times smaller than the 

acceleration of the Sun toward the center of the Galaxy, a ga i/c ~ a^/cx (di/8 kpc), where 

a>gai,o/ c — 6.5 x 10~ 19 s -1 . The probability map corrected for proper motions and for relative 

Galactic accelerations is shown in Figure 3. 

In contrast to the probability map for the simulated accelerated dataset (Figure 2), the 
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probability map based on the real data (Figure 3) does not show strong correlations in any 
direction and has on average much higher probabilities II of the null hypothesis. The lowest 
value of II in any direction in Figure 3 is 0.213, which is consistent with the null hypothesis. 
Therefore, there is no obvious evidence for a non-zero acceleration of the solar system relative 
to the MSP population. 

To establish the significance of this null result and to find the maximum acceleration 
allowed by the data on MSPs, we first created 1000 made-up datasets by using the positions 
of the real objects and randomly shuffling their P' new jP' values. For each direction on the 
sky (a, 5), these shuffled datasets provide a distribution of rank correlation coefficients of 
uncorrelated datasets. We then correct the dataset Pl )new /Pl for a putative acceleration 
a & directed toward (a, 5) and calculate the rank correlation coefficient for this corrected 
dataset. We keep increasing the value of a until the data start showing the presence of a 
correlation in a statistical sense, i.e., when the rank correlation coefficient reaches its top 
5% value as determined from the shuffled datasets (which are presumed uncorrelated with 
the position on the sky). When such value of a & is reached, we have overcorrected for the 
acceleration, and this value is ruled out by the data at the 95% level. We calculate these 
limiting values of acceleration for each direction on the sky and show them in Figure 4. 44% of 
the sky has limiting accelerations a Q /c < 10 -18 s -1 , 99% of the sky has limiting accelerations 
a©/c < 2 x 10 -18 s _1 and 100% of the sky has limiting accelerations a /c < 3 x 10 -18 s _1 . 

None of the conclusions of this Section change if we use the observed values of P'/P' 
without correcting for the proper motion effects and the relative Galactic accelerations. The 
corresponding probability map appears different in detail from the one shown in Figure 3, 
but it has a similar value of (II) and min(II) and shows no sign of a correlation. The limiting 
accelerations are also similar to the values found in the preceding paragraph. This finding 
is consistent with the fact that the limiting accelerations are at least a few times larger 
than the typical values of these corrections. This also means that the uncertainties in the 
distances to MSPs, however large, do not seriously affect the statistical method. The method 
is ultimately limited by the relatively broad distribution of the intrinsic values of P/P, as 
illustrated by equation (5). 

3. Individual objects 

Unlike the case of MSPs, where the intrinsic values of P/P are unknown, in some 
objects these values are known from theoretical arguments. Since the relative acceleration is 
constrained by the difference between the observed and the theoretical period change rates 
{Pobs — Pth)/P, both P obs and P t h must be measured or constrained with enough precision 
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to produce an interesting estimate of acceleration. Therefore, we consider only objects for 
which P obs /P < KT 17 s- 1 or so. 

The first group of such objects is short-period binary systems containing a pulsar (see 
Stairs 2004 for a recent review). From this list we excluded all objects found in globular 
clusters, which might be affected by peculiar accelerations associated with the cluster poten- 
tial. Of the remaining binary systems there are five in which the derivative P& of the orbital 
period Pb has been measured (Table 1), including the most recent measurement for the bi- 
nary pulsar PSR J0737-3039 (Burgay et al. 2003; Lyne et al. 2004; Kramer et al. 2005). In 
addition, there are five more binary systems for which interesting observational constraints 
exist on the orbital period decay. These ten systems are listed in Table 1 grouped by the 
accuracy of the constraint that they provide (the six objects on top of the Table form the 
higher accuracy 'group 1'). In these objects, the intrinsic value of Pb can be calculated under 
the assumptions that the orbital period decay is due to emission of gravitational waves and 
that this emission is correctly described by general relativity. 

The second group of objects with potentially known intrinsic P/P is pulsating white 
dwarfs. In this case P refers to the period of any of the oscillation modes. Pulsating 
white dwarfs come in three types, based on their effective temperature (Mukadam et al. 
2003); nearly a hundred objects are known (SIMBAD 2 ; Mukadam et al. 2004). While the 
typical periods of the oscillation modes are a few hundred seconds in all three classes, the 
period derivatives range from 10~ n to 10~ 15 , decreasing for cooler objects (Sullivan 1998). 
Because of our precision requirements, only the coolest pulsating white dwarfs (typically 
spectroscopically classified as hydrogen-rich, or DA) are suitable for acceleration constraints. 
In these stars the intrinsic value P is determined by the rate of cooling, while other factors 
such as gravitational contraction are unimportant (Kepler et al. 2000). Cooling of DA white 
dwarfs has been modeled leading to estimates of the intrinsic P/P (Bradley 1998). However, 
the period derivative has been measured only for a DA white dwarf G117— B15A, which is 
the only one that we include in our analysis (Kepler et al. 2000). For some other objects 
(e.g., the prototype of the class ZZ Ceti, Mukadam et al. 2003) upper limits on P b. s exist, 
but they are less restrictive than the actual measurement for G117— B15A. 

Kinematic corrections due to proper motions P pm ± a pm and due to relative Galactic 
accelerations P ga i ±(J ga i need to be taken into account for all objects that we use, in order to 
meaningfully compare the observed P Q b s ± <y bs and the theoretical Pth ± &th values. Timing 
studies listed above and in Table 1 typically include estimates of these corrections, unless the 
authors believe the corrections are unimportant. In one case (PSR B1913+16, Damour & 
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Taylor 1991) the kinematic correction also includes accelerations induced on the sun and the 
pulsar by their respective neighbor stars (with typical distances > 1 pc) and giant molecular 
clouds, but these authors find that these accelerations are tiny compared to the Galactic 
acceleration and the proper motion effect. 

We calculated the Galactic acceleration corrections using the Milky Way potential in the 
form suggested by P90. The errors of these corrections were calculated including both the 
uncertainty in the distance and the uncertainty of the model of the potential. The uncertainty 
in the model can be estimated by comparing the P90 potential to other suggested forms of 
potential, for example, the most recent update by Dehnen & Binney (1998). For each object, 
the relative error in the acceleration component parallel to the Galactic plane is calculated 
by adding in quadrature the relative error in the distance between the sun and this object 
and 15%. This last value comprises the uncertainties of the potential model in the plane 
(< 10%, Sun & Han 2004), the uncertainty in the solar Galactic velocity (< 10%) and 
the uncertainty in the distance to the Galactic center (< 10%). The relative error in the 
acceleration component perpendicular to the plane is obtained by adding in quadrature the 
relative distance error and 25% which is dominated by the uncertainty in the potential model 
(Sun & Han 2004). 

From equation (2), the acceleration of the Solar system in any direction can be con- 
strained from each object with a known P b s and P^: 

~P = Fobs ~ Pth ~ Ppm ~ Pgal ± + < + tfm + °gal = AP ± O toU (7) 

where the unit vector n points toward the object that we are using. If only one object is 
available, the acceleration is unconstrained in the plane perpendicular to the line of sight. 
Therefore, the constraint on the acceleration can be improved by combining data on several 
objects in different parts of the sky. 

In Table 1 we list all the values that enter equation (7) and their errors for the ten pulsars 
in binary systems and one DA white dwarf used in our analysis. The objects PSR B1913+16 
and PSR J1713+0747 provide the most interesting constraints on the acceleration. For 
PSR B1913+16, the accuracy improved by about an order of magnitude in the last 13 years 
(cf. Weisberg & Taylor 2004 and Damour & Taylor 1991). This object has been timed for 
30 years and has been used to constrain the acceleration of the solar system all by itself 
(Thornburg 1985). In the object PSR J1713+0747 the orbital period decay has not been 
detected, but the available upper limit on the orbital decay based on 12 years of observations 
combined with the long orbital period of this system (68 days) makes it useful for our analysis. 

The best-fit acceleration in the direction (a, 5) from object i is given by a Q /c = — (AP± 
&tot,i)/(Pi cos9i) where 6>j is the angle between (a, 5) and the direction to the object. There 
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is no evidence for a non-zero acceleration at the 2a level for any object. In Figure 5 we 
show the values of accelerations ruled out at a 2a level for each direction on the sky based 
on the observations of PSR B1913+16 alone. The least constrained directions are those 
perpendicular to the direction to the pulsar. The constraints in these regions are significantly 
improved when we include another five objects from the high-accuracy group, as shown in 
Figure 6. In this Figure, for each position on the sky we calculated the error- weighted mean 
and the error-weighted variance of the accelerations obtained using equation (7) for the six 
pulsars and plotted (mean + 2 x variance 1 / 2 ). Adding the remaining five objects from Table 
1 has almost no effect on the acceleration constraints. Table 2 and Figure 7 summarize what 
values of acceleration of the solar system can be ruled out by different methods. 

The method described here assumes that the peculiar accelerations of each object due 
to nearby stars and clouds are negligible compared to their Galactic accelerations relative to 
the sun. All objects used in our calculation are either in the disk of the Galaxy or somewhat 
above it, rather than in the bulge, and most are within 1 kpc, so the estimate of the peculiar 
acceleration of the sun (Section 2.2) is likely to apply to all these objects as well. If so, the 
peculiar acceleration is negligible compared to the uncertainties introduced by the distance 
determination. 



4. Constraints on the background of gravitational waves 

The universe may be filled with gravitational waves, produced either by individual 
sources such as supernovae or merging stars, or at very early stages of the expansion shortly 
after the Big Bang (Allen 1997). A passing wave affects the propagation of electromag- 
netic emission and produces a difference between the observed and the emitted frequencies. 
Specifically, the frequency change due to gravity waves is (Mashhoon 1979) 

^obs ^em / \ 
Vem - Vobs, (8) 

where v is the frequency and rj <C 1 is the perturbation of the metric due to gravitational 
waves at the position of the source (v em , r] em ) and of the observer {v ohs , r] obs ). Similarly, 
the observed rotational period and period change of a pulsar are different from the intrinsic 
values in the presence of gravitational waves (Sazhin 1978; Detweiler 1979). Of course, in 
order to use these differences to compute the effects of the gravitational waves, we now have 
to assume that there is no unknown contribution to the observed period change (i.e., from a 
massive companion). 

For example, gravitational waves with frequencies / > 1/T, where T = 10 — 30 years 
is the span of pulsar timing, introduce noise to the timing data. With some assumptions 
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about the spectral energy distribution of the waves, one can constrain the amplitude of 
the gravitational waves (and therefore their energy density) from observations of individual 
pulsars (Detweiler 1979; Romani & Taylor 1983; Kaspi et al. 1994; Thorsett & Dewey 1996; 
McHugh et al. 1996) or pulsar pairs using cross-correlation of timing noise from two objects 
(Hellings & Downs 1983). This method is sensitive to frequencies 1/T < / < 1/r, where r is 
a typical spacing between consecutive measurements, and the most stringent available upper 
limit on the energy density in such waves, expressed as a fraction of the critical density of 
the universe, is Vt GW h 2 < 1.0 x 1(T 8 (Thorsett & Dewey 1996; McHugh et al. 1996), where 
H = h x 100 km s" 1 Mpc -1 is the Hubble constant. 

In this paper we have used the values of the period and the period derivative which 
are the best-fit timing solutions over the entire span of available data, thereby discarding 
possible variations with frequencies / > 1/T. If we now focus on stochastic waves with 
frequencies c/d < / < 1/T (where d is the typical distance to the pulsars from the sun), 
then the r] em terms of equation (8) are uncorrelated for different pulsars, but the second term 
of this equation provides a common change of periods and period derivatives for all objects. 
The difference between the expected and the measured values of the period derivatives 
in individual objects, such as those listed in Table 1, can be directly used to constrain 
the energy density in gravitational waves in this frequency regime (Bertotti et al. 1983; 
Taylor & Weisberg 1989; Thorsett & Dewey 1996), giving Vt GW h 2 < 0.04 (Thorsett & Dewey 
1996). The accuracy of the AP/P measurement for PSR B1913+16 (which dominated this 
constraint) has improved by about a factor of three since the publication by Thorsett & 
Dewey (1996), so the current upper limit on Vl GW h 2 oc (AP/P) 2 is approximately an order 
of magnitude better. In contrast to the case of the peculiar solar acceleration, the effect of 
gravitational waves in this frequency range on timing measurements does not depend on the 
position on the sky, so it cannot be extracted statistically from the timing data as we did in 
Section 2 for the solar acceleration. 

Finally, for very low frequency gravitational waves (/ < c/d) the size of the 'detector' 
(the solar system and the pulsars) is smaller than the wavelength of the waves, and it can 
be shown that, to first order in 77, the change in the observed period of the signals is 

P'-P d ^ 

where rii are the spatial components of the unit vector in the direction to the pulsar and 
rjij are the spatial components of the metric perturbation taken at the observer's position. 
This regime can be probed both with individual objects (by taking a time derivative of 
eq. 9 and using limits on AP from Table 1) and in the statistical sense using MSPs since 
the effect is expected to increase with distance. In the latter case, one must be aware 
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that a systematic dependence on the distance may be introduced when applying the proper 
motion correction to the observed values P[j 'Pi, since /ij are typically not known for more 
distant objects. For the subset of MSPs with measured proper motions, the correlation of 
Pi/ Pi versus di is present at less than 70% confidence level (using the rank correlation test). 
For individual objects, the amplitude of the metric perturbation is best constrained with 
PSR J1713+0747 and PSR B1913+16 and is consistent with within l.lcr. Using equation 
(9) and the expression for energy density of gravitational waves (e.g., Detweiler 1979) we 
can place an upper limit on the energy density of very low frequency waves: 

With the best accuracy of Utot/ P — 1 x 10~ 19 s _1 for objects in Table 1, equation (10) implies 
ftlfh 2 < 0.004 (c/df) 2 (95% confidence), where f2/d/ is the density relative to the critical 
density in the frequency range (/, / + d/). 



5. Conclusions 

In the presence of an undiscovered massive companion, the known solar system barycen- 
ter would experience an acceleration which could be detected in observations of stable as- 
tronomical clocks. No such acceleration has been detected so far. The best limits on the 
acceleration are provided by the sources for which the intrinsic period change rate can be 
somehow estimated and then compared with the observed period change. Examples of such 
systems are pulsars in binary systems in which the orbital period decays due to emission 
of gravitational waves, and white dwarfs in which the pulsation periods change as the star 
cools. Alternatively, statistical limits on the acceleration can be provided by a group of 
objects of the same type distributed across the sky, even if the intrinsic period change rate 
of these objects is unknown. An example is millisecond pulsars, which provide sensitivity to 
the acceleration similar to that obtained using pulsars in binary systems and white dwarfs. 
The limits on the acceleration of the solar system barycenter that we obtained by different 
methods are summarized in Table 2. 

The acceleration exerted by the companion on the solar system is a Q /c = Gm p /(cd p ), 
where d p is the current distance to the companion and m p is its mass. Therefore, the upper 
limit on the mass of the companion as a function of the upper limit on the acceleration is 

m p = 0.2QMjX ( gg/£ - ) (— ^-Y. (11) 

p J \5 x 10" 19 s- V V 100 Au / 
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Companions with masses higher than this can be ruled out at a 2a confidence over about 
76% of the sky (Table 2). Our limit on the acceleration is about an order of magnitude better 
than the one available 15 years ago (Tremaine 1990) because of the improved accuracy of the 
timing observations of PSR B1913+16 and because several other objects became available. 

Timing of pulsars has been widely used to place constraints on the energy density of 
low-frequency gravitational waves. In Section 4 we summarized these methods and outlined 
the differences between the effects of a companion and of the gravitational waves on pulsar 
timing. 

The observational accuracy will improve with time for many objects, orbital period 
decay will likely be measured for many pulsars in binary systems, and more millisecond 
pulsars will be discovered. For some objects used in our analysis further observations are 
unlikely to sharpen the limit on the acceleration since the current constraint that they 
provide is determined mostly by the uncertainties in the distances to pulsars and in models 
of the potential of the Galaxy. This limit has now been reached for PSR B1913+16 and 
PSR B1534+12 (Table 1), and there is little prospect for more accurate independent distance 
measurements (Weisberg & Taylor 2004). The limit on the acceleration is significantly 
tightened by adding objects which are closer and have a parallax measurement with better 
than 10% precision, such as PSR J1737+0747 (the acceleration constraint provided by this 
object is still dominated by the timing precision and will therefore improve with time). 
The statistical method using MSPs is currently limited by the intrinsic distribution of the 
period derivatives and is unlikely to rapidly improve until much larger samples of MSPs 
become available. Accelerations exerted by field objects and the uncertainties in solar system 
ephemerides are negligible in both methods. 

Constraining the acceleration of the solar system using individual objects like pulsars 
in binary systems and white dwarfs relies on the theoretical arguments that are used to 
calculate the intrinsic period derivative. The reversed sequence has been used for some of 
these objects. For example, PSR B1913+16 was used to confirm the prediction of general 
relativity for the rate of emission of gravitational waves assuming that the solar acceleration 
is negligible (Damour k Taylor 1991; Weisberg & Taylor 2004). Similarly, G117-B15A has 
been used to constrain the variation of the gravitational constant G (Benvenuto et al. 2004; 
Biesiada & Malec 2004) and other exotic physics (Biesiada & Malec 2002). 

The acceleration of the solar system barycenter is not the only constraint on the mass 
and the distance of the companion; other methods have been described in detail by Tremaine 
(1990) and Hogg et al. (1991). One of the most sensitive of these is constraining the tidal 
force exerted by a distant companion (oc m p /dp) using the ephemerides of known planets 
and comets. Depending on the assumed accuracy of the measured ephemerides, the masses 
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that can be ruled out by the acceleration method and by the tidal method are similar for 
dp ~ 400 AU or so, but the acceleration method becomes more sensitive if the companion is 
at a larger distance (the limiting m p is oc d 2 v rather than d p ). Furthermore, the tidal method 
is much less sensitive away from the ecliptic plane, whereas the acceleration method can 
provide a uniform sensitivity over the entire sky On the downside, the acceleration method 
is not sensitive to an extended symmetrical disk such as a massive Kuiper belt, whereas the 
tidal method can detect such structures, though with reduced sensitivity. 

The Voyager and Pioneer spacecraft have provided distances to Uranus and Neptune 
with unprecedented precision using ranging measurements during flybys (Anderson et al. 
1989, 1995b). However, the high accuracy of these measurements does not translate di- 
rectly into a significant improvement in the ephemerides of the planets, since the ranging 
measurements constrain only the radial component of the position, whereas the transverse 
component is still determined using optical astro metric data. As a result, the constraints 
on the unseen mass in the solar system derived from the planet ephemerides that incorpo- 
rate ranging data (Anderson et al. 1995b) are not very different from those derived from 
the optical ephemerides alone (Hogg et al. 1991). In addition to providing planet distances 
and masses, spacecraft themselves serve as accurate gravity probes. The two Pioneers are 
especially suitable (Anderson et al. 2002) as they did not require many reorientation ma- 
neuvers (which complicate modeling of trajectories). An explicit search for Planet X in 
the Pioneer tracking data yielded a null result (Anderson 1988), with a typical sensitivity 
a/c< (0.7—1.4) x 10~ 18 s _1 (2a upper limit on the gravitational acceleration of the spacecraft 
relative to the expected value; Anderson et al. 1995a, 2002), about an order of magnitude 
less sensitive than the constraint based on the ephemerides of the outer planets. Recently 
Anderson et al. (2002) have argued that Pioneers exhibit a constant anomalous acceleration 
of a/c ~ (2.9 ± 0.4) x 10~ 18 s -1 , but this is probably due to non-gravitational effects. 

Yet another constraint on the mass and the position of the companion can be derived 
based on the fairly accurate alignment of the spin of the sun and the orbital angular momenta 
of planets implying that there is no massive companion on a highly inclined orbit that would 
have forced the planetary orbits out of this configuration (Goldreich & Ward 1972; Hogg et 
al. 1991). For a high-inclination companion at a large distance (d p > 100 AU), this constraint 
is considerably more sensitive than the acceleration method. 

A planetary companion to the solar system would reflect solar light and therefore could 
be detected in the optical or near-infrared by direct observations; furthermore, its own 
thermal radiation would be detectable in the mid- to far-infrared. Neither the optical surveys 
by Tombaugh (1961) and Kowal (1989), nor the IRAS All-Sky Survey (Neugebauer et al. 
1984) discovered any massive companions (other than Pluto). The optical surveys used 
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a baseline of a few days to detect large-parallax objects, resulting in a typical distance 
sensitivity d p < 400 AU, whereas the IRAS survey was sensitive out to 400 — 2500 AU, 
depending on the position on the sky. Neither survey covered 100% of the sky, and the flux 
and parallax sensitivities varied across the sky in all three cases, as described in detail by 
Kowal (1989) and Hogg et al. (1991). In what follows we use Tombaugh's and IRAS surveys, 
as they had the largest sky coverage. 

Modulo incompleteness, Tombaugh's and IRAS surveys had sufficient flux sensitivity to 
detect a rock-ice, low-albedo object with the mass given by equation (11) out to about d p = 90 
AU. For the calculation of the limiting distances for a rock-ice planet we followed Hogg et al. 
(1991) in assuming an albedo of 0.02, the limiting visual magnitude of Tombaugh's survey 
of 15.2, temperature of 280 K (d p /l AU)^ 1 / 2 (determined exclusively by solar heating), mass 
density of 2 g cm -3 and the limiting flux of IRAS at lOOyum of 1 Jy. If the unseen companion 
is similar to giant planets rather than a rocky body, then its thermal radiation is dominated 
by its evolutionary cooling rather than by solar heating, and its radius is almost independent 
of mass. We used cooling models of brown dwarfs by Allard et al. (2001) which cover the 
range of masses from 1 Mj to the hydrogen burning limit 0.08 M . In these models the 
temperature is a function of mass and age, and the isochrone at 5 Gyr can be approximated 
as 87 K (rrip/Mj) - 57 . A gas giant with mass given by equation (11) would be detectable by 
IRAS in the 60 and 100/xm bands farther than about 70 AU out to the distances when the 
parallax of such object would be undetectable. The optical emission from a giant planet is 
reflected sunlight, with typical albedo of up to 0.5 — 0.6, and is almost independent of its 
mass (insofar as the radius is independent of the mass). A giant planet with Jupiter's radius 
and albedo of 0.5 would be detectable in Tombaugh's survey out to 290 AU. 

To summarize the previous paragraphs, the acceleration method is more sensitive to 
distant (d p > 300 AU) gas giants or rock-ice bodies than Tombaugh's survey, but it is not as 
sensitive as IRAS survey to a gas giant at d p > 70 AU, out to about 400 AU, where IRAS 
parallax sensitivity becomes an issue. An alternative exotic possibility is that the companion 
is a dark compact object (a neutron star or a black hole) which could escape detectability 
in either survey. 

A mass given by equation (11) may be detectable in a future high-precision astrometric 
survey through microlensing of background stars (Cowling 1983), but would have been unde- 
tectable by Hipparcos. The minimum and the maximum masses detectable by this method 
are limited by the time coverage of the survey, since the microlensing event should be caught 
by at least one survey observation, but proceed faster than the total duration of the survey. 
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The resulting range of masses that can be detected at distances d p < 2500 AU is 

'"^(^O^^K^fe)' (12) 

where (3 is the astrometric sensitivity, r is the total duration of the mission and n ~ 70 is the 
typical number of survey scans of a given field. The numerical values of these parameters are 
given for the GAIA survey (to V — 15 mag). A more detailed discussion of the sensitivity 
of GAIA to distant companions is given by Gaudi & Bloom (2005). 
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Fig. 1. — Sky distribution of MSPs used for the analysis. 48 MSPs are marked with circles 
(filled for positive declinations and empty for negative declinations). For reference, the 
position of the Galactic center is indicated by a star. The same coordinate grid is used in 
all subsequent sky maps; the circles are lines of constant declinations, and right ascensions 
increase counter-clockwise. Equal solid angles are mapped into equal areas on this projection. 
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negative declinations 



0.25 0.5 0.75 1 

probability of the null hypothesis 



Fig. 2. — The probability map for a simulated population of pulsars with the same positions 
as observed MSPs, but with P = 0.007 s and P = 2.4 x 10~ 20 , if the solar system accelerates 
relative to this population with a Q /c = 10~ 18 s _1 in the direction a Q = 180, 5 Q = 0. Each 
area on the sky around the direction (a, 5) is shaded according to the probability that P[jP[ 
is not correlated with cos 9i relative to this direction. The darker the shading, the smaller the 
probability of the null hypothesis (no correlation), and the stronger is the correlation in this 
direction. The acceleration component is clearly visible (the light-colored band is the plane 
of the least correlation, which lies perpendicular to the direction of the acceleration); this 
highly significant detection of an acceleration has (II) = 0.113, min(II) ~ 0, max(II) ~ 1. 
The acceleration that can be retrieved based on the probability map is slightly offset relative 
to the input acceleration; this is due to the inhomogeneous distribution of objects on the 
sky. By construction, probability maps are invariant with respect to reversing the direction 
of the acceleration, n(a, 5) = U(ir + a, —5). The coordinate system and the legend are the 
same as in Figure 1. 




Fig. 3. — The probability map for the real data on MSPs corrected for the proper motion 
effect and the relative accelerations in the Galactic potential field. On this plot max(n) = 
0.996, min(n) = 0.213, (n) = 0.610. No acceleration signature is present in this map, as 
indicated by the high values of II (lighter shading compared to the previous Figure) and by 
the lack of reflection symmetry about a plane. The coordinate system and the legend are 
the same as in Figure 1. 
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Fig. 4. — Values of the acceleration that can be ruled out at a 95% confidence level in each 
direction based on timing of MSPs. The data have been corrected for proper motions and 
Galactic accelerations. The coordinate system and the legend are the same as in Figure 1. 
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Fig. 5. — Values of the acceleration of the solar system that can be ruled out at a la 
confidence level in each direction based on the observations of PSR B1913+16 alone. The 
most poorly constrained directions are those perpendicular to the line of sight to the pulsar, 
hence the symmetric appearance. The coordinate system and the legend are the same as in 
Figure 1. 
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Fig. 6. — Values of the acceleration that can be ruled out at a 2a confidence level in each 
direction based on the observations of six 'group 1' objects from Table 1. Adding the 
remaining five objects has almost no effect on the acceleration constrains. The coordinate 
system and the legend are the same as in Figure 1. 
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Fig. 7. — The values of acceleration that can be ruled out at 95% level by different methods: 
statistical analysis using MSPs (solid line), timing of PSR Bf9f3+16 alone (dotted line), 
timing of the six 'group 1' pulsars combined (dashed line). In the latter case, adding the 
remaining five objects in Table 1 has almost no effect on the acceleration constraints. This 
Figure is the graphical complement to Table 2. 



Table 1. Constraints on the acceleration from pulsars in binary systems and pulsating white dwarfs 



name 


a 


8 


Pobs/P 


crobs/P 


Pth/P 


Tth/P 


Ppm 1 ' P 


Opm/P Pgal/P 


a gal /P AP/P (Ttot/P 


jn, mas/yr 


d, pc 


I 


b 


comp. 


ref. 


PSR J0437-4715 


69.32 


-47.25 


7.3 


0.4 


-0.00073 


0.00007 


6.71 


0.14 


-0.065 


0.015 


0.69 


0.43 


140.892(6) 


139(3) 


253.39 


-41.96 


WD 


1 


PSR J1012+5307 


153.14 


53.12 





1.9 


-0.21 


0.05 


1.3 


0.2 


-0.11 


0.04 


-1.1 


1.9 


25.3(4) 


840(90) 


160.35 


50.86 


WD 


2 


PSR B1534+12 


234.29 


11.93 


-3.55 


0.36 


-5.28 


0.02 


1.10 


0.33 


-0.18 


0.07 


0.81 


0.49 


25.5(5) 


700(200) 


19.85 


48.34 


NS 


3 


PSR J1713+0747 


258.46 


7.79 





0.1 


-0.052 


-0.009 


0.11 


0.01 


-0.089 


0.032 


0.03 


0.11 


6.296(14) 


1100(100) 


28.75 


25.22 


WD 


4 


PSR B1855+09 


284.40 


9.72 


0.6 


1.2 


-0.00011 


0.00002 


0.085 


0.026 


0.0030 


0.0010 


0.5 


1.2 


6.16(11) 


900(300) 


42.29 


3.06 


WD 


5 


PSR B1913+16 


288.87 


16.11 


-86.66 


0.03 


-86.0867 


0.0007 


0.097 


0.027 


-0.54 


0.12 


-0.13 


0.12 


2.6(3) 


5900(940) 


49.97 


2.12 


NS 


6 


PSR J0621+1002 


95.34 


10.04 


0.0 


7.0 


-0.0011 


0.0003 


0.058 


0.017 


0.12 


0.04 


-0.2 


7.0 


3.5(6) 


1900(500) 200.57 


-2.01 


WD 


7 


PSR J0737-3039 


114.46 


-30.66 


-136 


9 


-141.21 


0.18 


<0.16 




-0.030 


0.010 


5.3 


9.1 


<10.5 


600(200) 


245.24 


-4.50 


2PSR 


8 


PSR J1141-6545 


175.28 


-65.76 


-25 


6 


-22 


6 


<1.3 


<0.7 


-0.3(?) 


0.05(?) 


-4 


9 


12(3) 


>3800 


295.79 


-3.86 


WD 


9 


PSR J2145-0750 


326.46 


-7.84 


0.0 


4.1 


-0.0009 


0.0002 


0.25 


0.07 


-0.15 


0.05 


-0.1 


4.1 


14.1(7) 


500(130) 


47.78 


-42.08 


WD 


10 


G117-B15A 


141.07 


35.28 


10.7 


6.5 


17.2 


4.6 


4.3 


1.6 


-0.04 


0.02 


-10.7 


8.2 


136(2) 


95(36) 


189.06 


45.46 


none 


11 



Note. — Right ascensions (a) and declinations (<$) are given in degrees for the epoch 2000; they are known with much better accuracy than given in the Table. All 
P/P and their uncertainties cr/P are in units 10 — 18 s — 1 . For pulsars in binary systems P and P refer to their orbital periods and orbital period derivatives and for 
the DA white dwarf G117— B15A these are the pulsation period and its derivative. AP = P bs — Pth ~ Ppm — Pgai, ^tot = °"o6s + a 1h + a pm + a g a V ^ or pulsars in 
binary systems, the error in the theoretical values is due to the error in the mass measurements, whereas for the white dwarf it is mostly due to the range of acceptable 
theoretical models, with mass being the primary source of uncertainty. Used in calculation of the relative Galactic acceleration and the proper motion correction are /i 
(proper motion in mas year -1 ), d (distance to the object in parsecs), I and b (Galactic longitude and latitude in degrees). Galactic acceleration corrections are calculated 
based on the P90 model for the potential of the Galaxy. The proper motion corrections and the Galactic acceleration corrections agree within the errors with those 
given in the literature, where available, except for the Galactic acceleration corrections for PSR J2145— 0750 (Lohmer et al. 2004) and PSR J1713+0747 (Splaver et al. 
2005) which are unimportant in both cases. The type of companion (WD - white dwarf, NS - neutron star, 2PSR - the system consists of two pulsars) is listed in the 
sccond-to-last column. Above the horizontal line are 'group 1' objects. 

References. — 1 - van Stratcn et al. (2001); 2 - Langc ct al. (2001); 3 - Stairs et al. (1998); 4 - Splaver et al. (2005); 5 - Kaspi et al. (1994); 6 - Damour & Taylor 
(1991); Weisberg & Taylor (2003, 2004); 7 - Splaver ct al. (2002); 8 - Lyne ct al. (2004); Kramer et al. (2005); 9 - Bailes et al. (2003); 10 - Lohmer et al. (2004); 11 - 
Bradley (1998); Kepler et al. (2000) 



Table 2. Constraints on the acceleration of the solar system by different methods 



method 


2 x 1CT 19 s- 1 


5 x icr 19 s- 1 


1 x 1CT 18 s- 1 


2 x icr 18 s- 1 


3 x 1(T 18 s- 1 


4 x 10~ 18 s- 1 


MSPs 


0% 


3% 


44% 


99% 


100% 


100% 


PSR B1913+16 


21% 


50% 


75% 


88% 


92% 


94% 


six pulsars in binary systems (group 1) 


25% 


76% 


94% 


99.5% 


100% 


100% 


all eleven objects from Table 1 


25% 


76% 


94% 


99.5% 


100% 


100% 



Note. — This table lists the fractions of the directions on the sky for which the acceleration of the solar system (given as 
Oq/c) can be ruled out at a given level with a 95% confidence (see also Figure 7). For example, using MSPs accelerations 
clq/c> 1 x 10~ 18 s _1 can be ruled out for 44% of the sky at 95% confidence. 



